"""
# -*- coding: utf-8 -*-
# @Time    : 2023/5/21 20:39
# @Author  : 王摇摆
# @FileName: SMO.py
# @Software: PyCharm
# @Blog    ：https://blog.csdn.net/weixin_44943389?type=blog
"""
import numpy as np

class SMO:
    """
    硬间隔支持向量机
    序列最小优化算法实现（Sequential minimal optimization/SMO）
    """

    def __init__(self, X, y):
        # 训练样本特征矩阵（N * p）
        self.X = X
        # 训练样本标签向量（N * EXP1）
        self.y = y
        # 拉格朗日乘子向量（N * EXP1）
        self.alpha = np.zeros(X.shape[0])
        # 误差向量，默认为负的标签向量（N * EXP1）
        self.errors = -y
        # 偏移量
        self.b = 0
        # 权重向量（p * EXP1）
        self.w = np.zeros(X.shape[1])
        # 代价值
        self.cost = -np.inf

    def fit(self, tol = 1e-6):
        """
        算法来自 John C. Platt 的论文
        https://www.microsoft.com/en-us/research/uploads/prod/1998/04/sequential-minimal-optimization.pdf
        """
        # 更新变化次数
        numChanged = 0
        # 是否检查全部
        examineAll = True
        while numChanged > 0 or examineAll:
            numChanged = 0
            if examineAll:
                for idx in range(X.shape[0]):
                    numChanged += self.update(idx)
            else:
                for idx in range(X.shape[0]):
                    if self.alpha[idx] <= 0:
                        continue
                    numChanged += self.update(idx)
            if examineAll:
                examineAll = False
            elif numChanged == 0:
                examineAll = True
            # 计算代价值
            cost = self.calcCost()
            if self.cost > 0:
                # 当代价值的变化小于容许的范围时结束算法
                rate = (cost - self.cost) / self.cost
                if rate <= tol:
                    break
            self.cost = cost

    def update(self, idx):
        """
        对下标为 idx 的拉格朗日乘子进行更新
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        # 检查当前拉格朗日乘子是否满足KKT条件，满足条件则直接返回 0
        if self.checkKKT(idx):
            return 0
        if len(alpha[(alpha != 0)]) > 1:
            # 按照｜E1 - E2｜最大的原则寻找第二个待优化的拉格朗日乘子的下标
            jdx = self.selectJdx(idx)
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 EXP1
            if self.updateAlpha(idx, jdx):
                return 1
        # 当未更新成功时，遍历不为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] == 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 EXP1
            if self.updateAlpha(idx, jdx):
                return 1
        # 当依然没有没有更新成功时，遍历为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] != 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 EXP1
            if self.updateAlpha(idx, jdx):
                return 1
        # 依然没有更新时返回 0
        return 0

    def selectJdx(self, idx):
        """
        寻找第二个待优化的拉格朗日乘子的下标
        """
        errors = self.errors
        if errors[idx] > 0:
            # 当误差项大于零时，选择误差向量中最小值的下标
            return np.argmin(errors)
        elif errors[idx] < 0:
            # 当误差项小于零时，选择误差向量中最大值的下标
            return np.argmax(errors)
        else:
            # 当误差项等于零时，选择误差向量中最大值和最小值的绝对值最大的下标
            minJdx = np.argmin(errors)
            maxJdx = np.argmax(errors)
            if max(np.abs(errors[minJdx]), np.abs(errors[maxJdx])) - errors[minJdx]:
                return minJdx
            else:
                return maxJdx

    def calcB(self):
        """
        计算偏移量
        分别计算每一个拉格朗日乘子不为零对应的偏移量后取其平均值
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_gt = alpha[alpha > 0]
        # 拉格朗日乘子向量中不为零的数量
        alpha_gt_len = len(alpha_gt)
        # 全部为零时直接返回 0
        if alpha_gt_len == 0:
            return 0
        # b = y - Wx，具体算法请参考文章中的说明
        X_gt = X[alpha > 0]
        y_gt = y[alpha > 0]
        alpha_gt_y = np.array(np.multiply(alpha_gt, y_gt)).reshape(-1, 1)
        s = np.sum(np.multiply(alpha_gt_y, X_gt), axis=0)
        return np.sum(y_gt - X_gt.dot(s)) / alpha_gt_len

    def calcCost(self):
        """
        计算代价值
        按照文章中的算法计算即可
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        cost = 0
        for idx in range(X.shape[0]):
            for jdx in range(X.shape[0]):
                cost = cost + (y[idx] * y[jdx] * X[idx].dot(X[jdx]) * alpha[idx] * alpha[jdx])
        return np.sum(alpha) - cost / 2

    def checkKKT(self, idx):
        """
        检查下标为 idx 的拉格朗日乘子是否满足 KKT 条件
        EXP1. alpha >= 0
        EXP2. y * f(x) - EXP1 >= 0
        EXP3. alpha * (y * f(x) - EXP1) = 0
        """
        y = self.y
        errors = self.errors
        alpha = self.alpha
        r = errors[idx] * y[idx]
        if (alpha[idx] > 0 and r == 0) or (alpha[idx] == 0 and r >= 0):
            return True
        return False

    def calcE(self):
        """
        计算误差向量
        E = f(x) - y
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_y = np.array(np.multiply(alpha, y)).reshape(-1, 1)
        errors = X.dot(X.T).dot(alpha_y).T[0] + self.b - y
        return errors

    def calcU(self, idx, jdx):
        """
        计算拉格朗日乘子上界，使两个待优化的拉格朗日乘子同时大于等于0
        按照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return 0
        else:
            return max(0.0, alpha[jdx] - alpha[idx])

    def calcV(self, idx, jdx):
        """
        计算拉格朗日乘子下界，使两个待优化的拉格朗日乘子同时大于等于0
        按照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return alpha[jdx] + alpha[idx]
        else:
            return np.inf

    def updateAlpha(self, idx, jdx):
        """
        对下标为 idx、jdx 的拉格朗日乘子进行更新
        按照文章中的算法计算即可
        """
        if idx == jdx:
            return False
        X = self.X
        y = self.y
        alpha = self.alpha
        errors = self.errors
        # idx 的误差项
        Ei = errors[idx]
        # jdx 的误差项
        Ej = errors[jdx]
        Kii = X[idx].dot(X[idx])
        Kjj = X[jdx].dot(X[jdx])
        Kij = X[idx].dot(X[jdx])
        # 计算 K
        K = Kii + Kjj - 2 * Kij
        oldAlphaIdx = alpha[idx]
        oldAlphaJdx = alpha[jdx]
        # 计算 jdx 的新拉格朗日乘子
        newAlphaJdx = oldAlphaJdx + y[jdx] * (Ei - Ej) / K
        U = self.calcU(idx, jdx)
        V = self.calcV(idx, jdx)
        if U == V:
            return False
        if newAlphaJdx < U:
            # 当新值超过上界时，修改其为上界
            newAlphaJdx = U
        if newAlphaJdx > V:
            # 当新值低于下界时，修改其为下界
            newAlphaJdx = V
        if oldAlphaJdx == newAlphaJdx:
            # 当新值与旧值相等时，判断为未更新，直接返回
            return False
        # 计算 idx 的新拉格朗日乘子
        newAlphaIdx = oldAlphaIdx + y[idx] * y[jdx] * (oldAlphaJdx - newAlphaJdx)
        # 更新权重向量
        self.w = self.w + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx] + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx]
        # 更新拉格朗日乘子向量
        alpha[idx] = newAlphaIdx
        alpha[jdx] = newAlphaJdx
        oldB = self.b
        # 重新计算偏移量
        self.b = self.calcB()
        # 重新计算误差向量
        eps = np.finfo(np.float32).eps
        newErrors = []
        for i in range(X.shape[0]):
            newError = errors[i] + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx].dot(X[i]) + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx].dot(X[i]) - oldB + self.b
            if np.abs(newError) < eps:
                newError = 0
            newErrors.append(newError)
        self.errors = newErrors
        return True



